This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.


# Load libraries
library(readr)
library(ggplot2)
library(nlme)
library(plyr)
library(gplots)
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
library(gstat)
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
library(sp)
require(MuMIn)
## Loading required package: MuMIn
# SE function
se<-function(x) sqrt(var(x, na.rm=T)/length(x))

ggplotRegression <- function (fit) {
  require(ggplot2)
  ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
    geom_point() +
    stat_smooth(method = "lm", col = "red") +
    labs(title = paste("Adj R2 = ",signif(summary(fit)$adj.r.squared, 5),
                       "Intercept =",signif(fit$coef[[1]],5 ),
                       " Slope =",signif(fit$coef[[2]], 5),
                       " P =",signif(summary(fit)$coef[2,4], 5)))
}
Sceloporus_specimens = read_delim("Sceloporus_specimens.txt", 
                                  "\t", escape_double = FALSE, trim_ws = TRUE)
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   SeriesID = col_character(),
##   Genus = col_character(),
##   Species = col_character(),
##   Date = col_character(),
##   Area = col_character(),
##   Site = col_character(),
##   Time = col_time(format = ""),
##   Sex = col_character(),
##   Age10g = col_character(),
##   VerbatimLocation = col_character()
## )
## See spec(...) for full column specifications.
# Select lizards >= 10 g only
Sceloporus_specimens_Age10g = Sceloporus_specimens[which(Sceloporus_specimens$Age10g=='Adult'),]
getwd()
## [1] "/Users/vanw/Documents/Projects/Sceloporus projects/Morph evn MS/Data archive"
# Descriptive stats
# SVL
summary(Sceloporus_specimens_Age10g$SVL_mm)  
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   64.00   71.75   76.00   76.66   81.00   93.00
# Mass_g
summary(Sceloporus_specimens_Age10g$Mass_g)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10.10   13.50   16.20   17.06   19.02   32.10
# Dorsal scale counts
summary(Sceloporus_specimens_Age10g$Scales_dorsal)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   36.00   39.50   41.00   41.46   44.00   48.00      57
# Dorsal color
summary(Sceloporus_specimens_Age10g$Color_Dorsal)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   18.98   34.78   37.75   37.52   41.21   48.78       7
# data summary function
data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}
# Table 1
df3.SVL <- data_summary(Sceloporus_specimens_Age10g, varname="SVL_mm",
                    groupnames=c("Site"))
write.table(df3.SVL, "df3.SVL.txt⁩", sep="\t")
df3.dorsal.scales <- data_summary(Sceloporus_specimens_Age10g, varname="Scales_dorsal", 
                                   groupnames=c("Site"))
write.table(df3.dorsal.scales, "df3.dorsal.scales.txt⁩", sep="\t")
df3.dorsal.color <- data_summary(Sceloporus_specimens_Age10g, varname="Color_Dorsal", 
                                  groupnames=c("Site"))
write.table(df3.dorsal.color, "df3.dorsal.color.txt⁩", sep="\t")
# SVL
par(mfrow=c(1,1))
hist(Sceloporus_specimens_Age10g$SVL_mm, breaks = 20)

# Summer
# Mean Temp of Warmest Quarter

par(mfrow=c(1,1))
par(mar=c(5,5,2,2))

library(lme4)
## Loading required package: Matrix
## Registered S3 method overwritten by 'lme4':
##   method         from 
##   predict.merMod MuMIn
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
reg1d.lmer = lmer(SVL_mm ~ BIO10_meanTwarmestQ + (1|BIO10_meanTwarmestQ), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1d.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO10_meanTwarmestQ + (1 | BIO10_meanTwarmestQ)
##    Data: Sceloporus_specimens_Age10g
## 
##      AIC      BIC   logLik deviance df.resid 
##   1073.6   1086.0   -532.8   1065.6      160 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0984 -0.5958 -0.1051  0.5966  2.5261 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  BIO10_meanTwarmestQ (Intercept) 17.47    4.180   
##  Residual                        30.91    5.559   
## Number of obs: 164, groups:  BIO10_meanTwarmestQ, 33
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)          74.4157     5.1266 31.5532   14.52 1.63e-15 ***
## BIO10_meanTwarmestQ   0.1550     0.2673 31.6970    0.58    0.566    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BIO10_mnTwQ -0.984
anova(reg1d.lmer)
plot(reg1d.lmer)

r.squaredGLMM(reg1d.lmer)
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help
## page.
##              R2m       R2c
## [1,] 0.005909285 0.3648981
library(ggplot2)
reg1d.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO10_meanTwarmestQ, y = SVL_mm)) + 
  geom_point()
reg1d.1 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter (°C)", y = "Snout-Vent Length (mm)")

# Max Temp Warmest Month

library(lme4)
library(lmerTest)
reg1e.lmer = lmer(SVL_mm ~ BIO5_maxTwarmestMo + (1|BIO5_maxTwarmestMo), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1e.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO5_maxTwarmestMo + (1 | BIO5_maxTwarmestMo)
##    Data: Sceloporus_specimens_Age10g
## 
##      AIC      BIC   logLik deviance df.resid 
##   1071.4   1083.8   -531.7   1063.4      160 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2225 -0.6141 -0.1145  0.6227  2.4730 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  BIO5_maxTwarmestMo (Intercept) 20.66    4.545   
##  Residual                       29.39    5.421   
## Number of obs: 164, groups:  BIO5_maxTwarmestMo, 36
## 
## Fixed effects:
##                    Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)         74.1033     6.6428 30.1156  11.155 3.21e-12 ***
## BIO5_maxTwarmestMo   0.1160     0.2228 30.5010   0.521    0.606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BIO5_mxTwrM -0.990
anova(reg1e.lmer)
plot(reg1e.lmer)

r.squaredGLMM(reg1e.lmer)
##              R2m       R2c
## [1,] 0.004640269 0.4155082
library(ggplot2)
reg1d.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO5_maxTwarmestMo, y = SVL_mm)) + 
  geom_point()
reg1d.1 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month (°C)", y = "Snout-Vent Length (mm)")

# Mean Temp of Coldest Quarter

library(lme4)
library(lmerTest)
reg1b.lmer = lmer(SVL_mm ~ BIO11_meanTcoldestQ + (1|BIO11_meanTcoldestQ), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1b.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO11_meanTcoldestQ + (1 | BIO11_meanTcoldestQ)
##    Data: Sceloporus_specimens_Age10g
## 
##      AIC      BIC   logLik deviance df.resid 
##   1073.9   1086.3   -532.9   1065.9      160 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.46586 -0.63245 -0.03986  0.62416  3.02117 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  BIO11_meanTcoldestQ (Intercept) 17.38    4.168   
##  Residual                        31.27    5.592   
## Number of obs: 164, groups:  BIO11_meanTcoldestQ, 33
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)          78.3629     1.2454 29.2064  62.921   <2e-16 ***
## BIO11_meanTcoldestQ  -0.4895     0.3630 28.9627  -1.348    0.188    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BIO11_mnTcQ -0.655
anova(reg1b.lmer)
plot(reg1b.lmer)

r.squaredGLMM(reg1b.lmer)
##             R2m       R2c
## [1,] 0.03352773 0.3787038
library(ggplot2)
reg1d.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO11_meanTcoldestQ, y = SVL_mm)) + 
  geom_point()
reg1d.1 + theme_classic() + labs(x = "Mean Temperature of Coldest Quarter (°C)", y = "Snout-Vent Length (mm)")

# Min temp coldest month

library(lme4)
library(lmerTest)
reg1c.lmer = lmer(SVL_mm ~ BIO6_minTcoldestMo + (1|BIO6_minTcoldestMo), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1c.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: SVL_mm ~ BIO6_minTcoldestMo + (1 | BIO6_minTcoldestMo)
##    Data: Sceloporus_specimens_Age10g
## 
##      AIC      BIC   logLik deviance df.resid 
##   1066.3   1078.7   -529.2   1058.3      160 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.92480 -0.63242 -0.08504  0.60222  2.45331 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  BIO6_minTcoldestMo (Intercept) 13.48    3.671   
##  Residual                       30.49    5.521   
## Number of obs: 164, groups:  BIO6_minTcoldestMo, 33
## 
## Fixed effects:
##                    Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)         73.1123     1.8338 29.2611   39.87   <2e-16 ***
## BIO6_minTcoldestMo  -0.8527     0.3452 28.7928   -2.47   0.0197 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BIO6_mnTclM 0.887
anova(reg1c.lmer)
plot(reg1c.lmer)

r.squaredGLMM(reg1c.lmer)
##             R2m       R2c
## [1,] 0.09776281 0.3743755
library(ggplot2)
reg1c.1 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO6_minTcoldestMo, y = SVL_mm)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "black")
reg1c.1 + theme_classic() + labs(x = "Minimum Temperature of Coldest Month (°C)", y = "Snout-Vent Length (mm)")

reg1c <- lm(Sceloporus_specimens_Age10g$SVL_mm~Sceloporus_specimens_Age10g$BIO6_minTcoldestMo)
summary(reg1c)
## 
## Call:
## lm(formula = Sceloporus_specimens_Age10g$SVL_mm ~ Sceloporus_specimens_Age10g$BIO6_minTcoldestMo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.8321  -5.3122  -0.2677   4.7100  15.8614 
## 
## Coefficients:
##                                                Estimate Std. Error t value
## (Intercept)                                     72.6509     1.1003  66.030
## Sceloporus_specimens_Age10g$BIO6_minTcoldestMo  -0.8467     0.2044  -4.143
##                                                Pr(>|t|)    
## (Intercept)                                     < 2e-16 ***
## Sceloporus_specimens_Age10g$BIO6_minTcoldestMo 5.51e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.679 on 162 degrees of freedom
## Multiple R-squared:  0.0958, Adjusted R-squared:  0.09022 
## F-statistic: 17.16 on 1 and 162 DF,  p-value: 5.507e-05

Spatially-explicit models

note: some values my differ slightly from publication values due to variation in random jitter function of gps points. seed set for all spatially explict models in this example for reproducibility.

# Load dataset
library(readr)
library(ggplot2)
library(nlme)
library(plyr)
library(gplots)

# SE function
se<-function(x) sqrt(var(x, na.rm=T)/length(x))
Sceloporus_VertNet = read_delim("Sceloporus_occidentalis_VertNet.txt", 
                                  "\t", escape_double = FALSE, trim_ws = TRUE)
## Parsed with column specification:
## cols(
##   SVL_mm = col_double(),
##   decimallatitude = col_double(),
##   decimallongitude = col_double(),
##   SeriesID = col_character(),
##   bio_5 = col_double(),
##   bio_6 = col_double(),
##   bio_10 = col_double(),
##   bio_11 = col_double()
## )
summary(Sceloporus_VertNet)
##      SVL_mm      decimallatitude decimallongitude   SeriesID        
##  Min.   : 24.0   Min.   :32.89   Min.   :-122.7   Length:752        
##  1st Qu.: 64.0   1st Qu.:37.35   1st Qu.:-119.7   Class :character  
##  Median : 73.0   Median :37.75   Median :-119.5   Mode  :character  
##  Mean   : 70.6   Mean   :37.65   Mean   :-119.3                     
##  3rd Qu.: 79.0   3rd Qu.:37.93   3rd Qu.:-118.6                     
##  Max.   :109.0   Max.   :47.59   Max.   :-115.4                     
##      bio_5           bio_6             bio_10          bio_11      
##  Min.   :198.0   Min.   :-104.00   Min.   :110.0   Min.   :-29.00  
##  1st Qu.:241.0   1st Qu.: -73.00   1st Qu.:144.0   1st Qu.: -9.00  
##  Median :285.5   Median : -54.00   Median :181.0   Median : 15.00  
##  Mean   :277.5   Mean   : -47.67   Mean   :175.7   Mean   : 17.71  
##  3rd Qu.:307.0   3rd Qu.: -23.00   3rd Qu.:200.0   3rd Qu.: 37.00  
##  Max.   :364.0   Max.   :  58.00   Max.   :260.0   Max.   :122.00
# Adjust temperature scale
Sceloporus_VertNet$bio_5 = Sceloporus_VertNet$bio_5/10
Sceloporus_VertNet$bio_6 = Sceloporus_VertNet$bio_6/10
Sceloporus_VertNet$bio_10 = Sceloporus_VertNet$bio_10/10
Sceloporus_VertNet$bio_11 = Sceloporus_VertNet$bio_11/10
# Select lizards 64 g and larger only
Sceloporus_VertNet_65mm = Sceloporus_VertNet[which(Sceloporus_VertNet$SVL_mm>63.5),]
summary(Sceloporus_VertNet_65mm)
##      SVL_mm       decimallatitude decimallongitude   SeriesID        
##  Min.   : 64.00   Min.   :33.68   Min.   :-122.4   Length:577        
##  1st Qu.: 71.00   1st Qu.:37.32   1st Qu.:-119.7   Class :character  
##  Median : 76.00   Median :37.74   Median :-119.5   Mode  :character  
##  Mean   : 76.47   Mean   :37.64   Mean   :-119.2                     
##  3rd Qu.: 81.00   3rd Qu.:37.93   3rd Qu.:-118.6                     
##  Max.   :109.00   Max.   :47.59   Max.   :-115.4                     
##      bio_5           bio_6             bio_10          bio_11      
##  Min.   :19.80   Min.   :-10.400   Min.   :11.00   Min.   :-2.900  
##  1st Qu.:23.70   1st Qu.: -7.300   1st Qu.:14.30   1st Qu.:-1.000  
##  Median :27.80   Median : -5.800   Median :17.50   Median : 1.200  
##  Mean   :27.66   Mean   : -5.035   Mean   :17.46   Mean   : 1.555  
##  3rd Qu.:31.00   3rd Qu.: -2.800   3rd Qu.:20.00   3rd Qu.: 3.600  
##  Max.   :36.40   Max.   :  5.800   Max.   :26.00   Max.   :12.200
hist(Sceloporus_VertNet_65mm$SVL_mm)

Mean Temperature of Warmest Quarter (bio10)

library(rgdal)
## rgdal: version: 1.4-6, (SVN revision 841)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000
library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))

library(geoR)
## Warning in system2("/usr/bin/otool", c("-L", shQuote(DSO)), stdout = TRUE):
## running command ''/usr/bin/otool' -L '/Library/Frameworks/R.framework/
## Resources/library/tcltk/libs//tcltk.so'' had status 1
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
par(mfrow=c(1,1))
plot(v1)

library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
                            resamp=10, quiet = FALSE)
## 1  of  10 
2  of  10 
3  of  10 
4  of  10 
5  of  10 
6  of  10 
7  of  10 
8  of  10 
9  of  10 
10  of  10 
plot(ncf.scor)

mod.nocorr = lm(SVL_mm ~ bio_10, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
plot(v1)

# jitter GPS points 
set.seed(9999)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)
## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_10, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corExp(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corExp_1a = SVL_mm.corExp
summary(SVL_mm.corExp_1a)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3733.249 3755.021 -1861.624
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 52.0663289  0.5255908 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 80.68188  3.406257 23.686374  0.0000
## bio_10      -0.24457  0.154141 -1.586665  0.1131
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.92 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66309328 -0.68428635 -0.07113418  0.56200949  4.28756934 
## 
## Residual standard error: 7.751242 
## Degrees of freedom: 577 total; 575 residual
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_1a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_10, 
                     data = Sceloporus_VertNet_65mm, 
                     correlation = corGaus(form = ~j.easting + j.northing, 
                                           nugget = TRUE, 
                                           value = c(50, 0.5)))
SVL_mm.corGaus_1b = SVL_mm.corGaus
summary(SVL_mm.corGaus_1b)      
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3738.805 3760.577 -1864.402
## 
## Correlation Structure: Gaussian spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 31.4080885  0.5711817 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 81.69492 2.9314082 27.868831  0.0000
## bio_10      -0.29681 0.1365493 -2.173631  0.0301
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.937
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.72288021 -0.69083161 -0.09976393  0.57325596  4.37452618 
## 
## Residual standard error: 7.605604 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_1b, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corGaus.var)

# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_10, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corLin(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corLin_1c = SVL_mm.corLin
summary(SVL_mm.corLin_1c)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC    logLik
##   3734.05 3755.822 -1862.025
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.5094906  0.5646552 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 79.58240  3.348365 23.767543  0.0000
## bio_10      -0.18927  0.161908 -1.169009  0.2429
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.963
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.70130731 -0.70563209 -0.04864043  0.61073152  4.46851351 
## 
## Residual standard error: 7.434682 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_1c, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corLin.var)

# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_10, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corRatio(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corRatio_1d = SVL_mm.corRatio
summary(SVL_mm.corRatio_1d)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3738.134 3759.906 -1864.067
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 32.7784921  0.5579583 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 81.58379  3.091487 26.389823  0.0000
## bio_10      -0.29214  0.138558 -2.108464  0.0354
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.909
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.69715109 -0.68172480 -0.09456885  0.56885085  4.32730113 
## 
## Residual standard error: 7.692627 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_1d, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corRatio.var)

# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_10, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corSpher(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corSpher_1e = SVL_mm.corSpher
summary(SVL_mm.corSpher_1e)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3735.045 3756.817 -1862.523
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 90.8478140  0.5282386 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 81.36789  3.276044 24.837239  0.0000
## bio_10      -0.27635  0.152577 -1.811218  0.0706
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.938
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6882349 -0.6930102 -0.1013650  0.5508959  4.2819509 
## 
## Residual standard error: 7.75038 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_1e, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corSpher.var)

# calculate AIC to select best model
AIC(SVL_mm.corExp_1a, SVL_mm.corGaus_1b, SVL_mm.corLin_1c, SVL_mm.corRatio_1d, SVL_mm.corSpher_1e)
# SVL_mm.corExp is best
# SVL_mm.corSpher & SVL_mm.corLin are close
# best model
summary(SVL_mm.corExp_1a)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3733.249 3755.021 -1861.624
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 52.0663289  0.5255908 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 80.68188  3.406257 23.686374  0.0000
## bio_10      -0.24457  0.154141 -1.586665  0.1131
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.92 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66309328 -0.68428635 -0.07113418  0.56200949  4.28756934 
## 
## Residual standard error: 7.751242 
## Degrees of freedom: 577 total; 575 residual
# other good models
summary(SVL_mm.corLin_1c)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC    logLik
##   3734.05 3755.822 -1862.025
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.5094906  0.5646552 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 79.58240  3.348365 23.767543  0.0000
## bio_10      -0.18927  0.161908 -1.169009  0.2429
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.963
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.70130731 -0.70563209 -0.04864043  0.61073152  4.46851351 
## 
## Residual standard error: 7.434682 
## Degrees of freedom: 577 total; 575 residual
summary(SVL_mm.corSpher_1e)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_10 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3735.045 3756.817 -1862.523
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 90.8478140  0.5282386 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 81.36789  3.276044 24.837239  0.0000
## bio_10      -0.27635  0.152577 -1.811218  0.0706
## 
##  Correlation: 
##        (Intr)
## bio_10 -0.938
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6882349 -0.6930102 -0.1013650  0.5508959  4.2819509 
## 
## Residual standard error: 7.75038 
## Degrees of freedom: 577 total; 575 residual
# all give same interpretation
library(ggplot2)
plot3 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_10, y = SVL_mm)) + 
  geom_point()
plot3 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter", y = "Snout-Vent Length (mm)")

reg3 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_10)
summary(reg3)
## 
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_10)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.983  -5.341  -0.544   4.345  33.770 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     84.6603     1.4645  57.809  < 2e-16 ***
## Sceloporus_VertNet_65mm$bio_10  -0.4692     0.0821  -5.714 1.77e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.211 on 575 degrees of freedom
## Multiple R-squared:  0.05374,    Adjusted R-squared:  0.05209 
## F-statistic: 32.66 on 1 and 575 DF,  p-value: 1.768e-08
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_1a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

#plot
library(ggplot2)
plot3.1 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_10, y = SVL_mm)) + 
  geom_point()
plot3.1 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter", y = "Snout-Vent Length (mm)")

Max Temp of Warmest Month (bio_5)

library(rgdal)
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000
library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))

library(geoR)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
par(mfrow=c(1,1))
plot(v1)

library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
                            resamp=10, quiet = FALSE)
## 1  of  10 
2  of  10 
3  of  10 
4  of  10 
5  of  10 
6  of  10 
7  of  10 
8  of  10 
9  of  10 
10  of  10 
plot(ncf.scor)

mod.nocorr = lm(SVL_mm ~ bio_5, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
plot(v1)

# jitter GPS points 
set.seed(8888)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)
## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_5, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corExp(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corExp_2a = SVL_mm.corExp
summary(SVL_mm.corExp_2a)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC   logLik
##   3734.34 3756.112 -1862.17
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 55.0865638  0.5194307 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 81.57497  4.241857 19.230957  0.0000
## bio_5       -0.19039  0.129958 -1.465044  0.1435
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.946
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6321856 -0.6834408 -0.0421755  0.5776624  4.2676357 
## 
## Residual standard error: 7.804842 
## Degrees of freedom: 577 total; 575 residual
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_2a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_5, 
                     data = Sceloporus_VertNet_65mm, 
                     correlation = corGaus(form = ~j.easting + j.northing, 
                                           nugget = TRUE, 
                                           value = c(50, 0.5)))
SVL_mm.corGaus_2b = SVL_mm.corGaus
summary(SVL_mm.corGaus_2b)      
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##        AIC     BIC    logLik
##   3740.048 3761.82 -1865.024
## 
## Correlation Structure: Gaussian spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 31.9129845  0.5650735 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 82.77308  3.697749 22.384722  0.0000
## bio_5       -0.22993  0.115757 -1.986289  0.0475
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.96 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.69063367 -0.69890215 -0.07982741  0.57979998  4.35722315 
## 
## Residual standard error: 7.649751 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_2b, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corGaus.var)

# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_5, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corLin(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corLin_2c = SVL_mm.corLin
summary(SVL_mm.corLin_2c)  
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3735.683 3757.455 -1862.841
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.7132778  0.5630984 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 79.40055  4.317961 18.388433  0.0000
## bio_5       -0.11801  0.138808 -0.850169  0.3956
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.978
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66430333 -0.70442457 -0.04417102  0.63742255  4.46110240 
## 
## Residual standard error: 7.452414 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_2c, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corLin.var)

# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_5, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corRatio(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corRatio_2d = SVL_mm.corRatio
summary(SVL_mm.corRatio_2d)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3739.165 3760.937 -1864.583
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.5808965  0.5514612 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 82.72663  3.833582 21.579463    0.00
## bio_5       -0.22947  0.116857 -1.963699    0.05
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.939
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66625538 -0.68658515 -0.07497444  0.57688858  4.31019769 
## 
## Residual standard error: 7.740728 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_2d, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corRatio.var)

# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_5, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corSpher(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corSpher_2e = SVL_mm.corSpher
summary(SVL_mm.corSpher_2e)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC    logLik
##   3736.25 3758.022 -1863.125
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 91.7092542  0.5251655 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 82.31749  4.146448 19.852530  0.0000
## bio_5       -0.21249  0.129907 -1.635743  0.1024
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.961
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66136760 -0.70130724 -0.07443573  0.56013969  4.27534825 
## 
## Residual standard error: 7.776811 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_2e, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corSpher.var)

# calculate AIC to select best model
AIC(SVL_mm.corExp_2a, SVL_mm.corGaus_2b, SVL_mm.corLin_2c, SVL_mm.corRatio_2d, SVL_mm.corSpher_2e)
# SVL_mm.corExp
# best model
summary(SVL_mm.corExp_2a)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC   logLik
##   3734.34 3756.112 -1862.17
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 55.0865638  0.5194307 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 81.57497  4.241857 19.230957  0.0000
## bio_5       -0.19039  0.129958 -1.465044  0.1435
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.946
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6321856 -0.6834408 -0.0421755  0.5776624  4.2676357 
## 
## Residual standard error: 7.804842 
## Degrees of freedom: 577 total; 575 residual
# other models with support
summary(SVL_mm.corLin_2c)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3735.683 3757.455 -1862.841
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.7132778  0.5630984 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 79.40055  4.317961 18.388433  0.0000
## bio_5       -0.11801  0.138808 -0.850169  0.3956
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.978
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66430333 -0.70442457 -0.04417102  0.63742255  4.46110240 
## 
## Residual standard error: 7.452414 
## Degrees of freedom: 577 total; 575 residual
summary(SVL_mm.corSpher_2e)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_5 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC    logLik
##   3736.25 3758.022 -1863.125
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 91.7092542  0.5251655 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 82.31749  4.146448 19.852530  0.0000
## bio_5       -0.21249  0.129907 -1.635743  0.1024
## 
##  Correlation: 
##       (Intr)
## bio_5 -0.961
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66136760 -0.70130724 -0.07443573  0.56013969  4.27534825 
## 
## Residual standard error: 7.776811 
## Degrees of freedom: 577 total; 575 residual
# same interpretation
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_2a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

library(ggplot2)
plot4 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_5, y = SVL_mm)) + 
  geom_point()
plot4 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month", y = "Snout-Vent Length (mm)")

reg4 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_5)
summary(reg4)
## 
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_5)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.715  -5.440  -0.655   4.576  33.709 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   86.54793    1.97805  43.754  < 2e-16 ***
## Sceloporus_VertNet_65mm$bio_5 -0.36432    0.07066  -5.156 3.48e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.247 on 575 degrees of freedom
## Multiple R-squared:  0.04418,    Adjusted R-squared:  0.04252 
## F-statistic: 26.58 on 1 and 575 DF,  p-value: 3.484e-07

Winter

Mean Temp of Coldest Quarter (bio_11)

library(rgdal)
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000
library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))

library(geoR)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
par(mfrow=c(1,1))
plot(v1)

library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
                            resamp=10, quiet = FALSE)
## 1  of  10 
2  of  10 
3  of  10 
4  of  10 
5  of  10 
6  of  10 
7  of  10 
8  of  10 
9  of  10 
10  of  10 
plot(ncf.scor)

mod.nocorr = lm(SVL_mm ~ bio_11, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
plot(v1)

# jitter GPS points 
set.seed(9999)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)
## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_11, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corExp(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corExp_3a = SVL_mm.corExp
summary(SVL_mm.corExp_3a)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3729.706 3751.477 -1859.853
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 39.5944966  0.5678399 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.36239 1.3056903 59.25019  0.0000
## bio_11      -0.44154 0.1803243 -2.44860  0.0146
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.498
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.92580993 -0.72430551 -0.09372455  0.56376069  4.45828937 
## 
## Residual standard error: 7.442991 
## Degrees of freedom: 577 total; 575 residual
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_3a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_11, 
                     data = Sceloporus_VertNet_65mm, 
                     correlation = corGaus(form = ~j.easting + j.northing, 
                                           nugget = TRUE, 
                                           value = c(50, 0.5)))
SVL_mm.corGaus_3b = SVL_mm.corGaus
summary(SVL_mm.corGaus_3b)      
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3734.095 3755.867 -1862.048
## 
## Correlation Structure: Gaussian spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 29.3955433  0.6049416 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.50703 1.1169734 69.39022  0.0000
## bio_11      -0.48784 0.1601221 -3.04665  0.0024
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.516
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.97514173 -0.75594120 -0.09841634  0.55577894  4.49754792 
## 
## Residual standard error: 7.381887 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_3b, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corGaus.var)

# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_11, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corLin(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corLin_3c = SVL_mm.corLin
summary(SVL_mm.corLin_3c)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3728.694 3750.466 -1859.347
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.4939287  0.5926691 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.41861 1.0533765 73.49567  0.0000
## bio_11      -0.45026 0.1707887 -2.63636  0.0086
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.572
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.98416021 -0.74744628 -0.09743965  0.56886085  4.56579639 
## 
## Residual standard error: 7.262108 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_3c, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corLin.var)

# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_11, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corRatio(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corRatio_3d = SVL_mm.corRatio
summary(SVL_mm.corRatio_3d)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##       AIC      BIC    logLik
##   3734.23 3756.001 -1862.115
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##     range    nugget 
## 28.774088  0.594429 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.40568 1.3078337 59.18618  0.0000
## bio_11      -0.47775 0.1678583 -2.84615  0.0046
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.459
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.94268900 -0.73327320 -0.08573498  0.56136307  4.47032918 
## 
## Residual standard error: 7.441609 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_3d, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corRatio.var)

# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_11, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corSpher(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corSpher_3e = SVL_mm.corSpher
summary(SVL_mm.corSpher_3e)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3730.683 3752.455 -1860.342
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 53.5573820  0.5938406 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.46890 1.0923423 70.91998   0.000
## bio_11      -0.46464 0.1715471 -2.70851   0.007
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.562
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.99123398 -0.75453630 -0.09968822  0.56182758  4.55617772 
## 
## Residual standard error: 7.277445 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_3e, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corSpher.var)

# calculate AIC to select best model
AIC(SVL_mm.corExp_3a, SVL_mm.corGaus_3b, SVL_mm.corLin_3c, SVL_mm.corRatio_3d, SVL_mm.corSpher_3e)
#SVL_mm.corLin
summary(SVL_mm.corLin_3c)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3728.694 3750.466 -1859.347
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.4939287  0.5926691 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.41861 1.0533765 73.49567  0.0000
## bio_11      -0.45026 0.1707887 -2.63636  0.0086
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.572
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.98416021 -0.74744628 -0.09743965  0.56886085  4.56579639 
## 
## Residual standard error: 7.262108 
## Degrees of freedom: 577 total; 575 residual
#best support
# also look at these with approx. equal support
summary(SVL_mm.corExp_3a)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3729.706 3751.477 -1859.853
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 39.5944966  0.5678399 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.36239 1.3056903 59.25019  0.0000
## bio_11      -0.44154 0.1803243 -2.44860  0.0146
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.498
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.92580993 -0.72430551 -0.09372455  0.56376069  4.45828937 
## 
## Residual standard error: 7.442991 
## Degrees of freedom: 577 total; 575 residual
summary(SVL_mm.corSpher_3e)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_11 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3730.683 3752.455 -1860.342
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 53.5573820  0.5938406 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 77.46890 1.0923423 70.91998   0.000
## bio_11      -0.46464 0.1715471 -2.70851   0.007
## 
##  Correlation: 
##        (Intr)
## bio_11 -0.562
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.99123398 -0.75453630 -0.09968822  0.56182758  4.55617772 
## 
## Residual standard error: 7.277445 
## Degrees of freedom: 577 total; 575 residual
# same result
reg1 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_11)
summary(reg1)
## 
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_11)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.273  -5.166  -0.433   4.237  33.983 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    77.63003    0.32863 236.225  < 2e-16 ***
## Sceloporus_VertNet_65mm$bio_11 -0.74662    0.09529  -7.835 2.28e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.046 on 575 degrees of freedom
## Multiple R-squared:  0.09646,    Adjusted R-squared:  0.09489 
## F-statistic: 61.39 on 1 and 575 DF,  p-value: 2.276e-14
# abline(reg1, col = "black")
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_3a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

#plot
library(ggplot2)
plot1 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_11, y = SVL_mm)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "black")
plot1 + theme_classic() + labs(x = "Mean Temperature of Coldest Quarter", y = "Snout-vent Length (mm)")

Min Temp of Coldest Month (bio_6)

library(rgdal)
xy = cbind(Sceloporus_VertNet_65mm$decimallongitude, Sceloporus_VertNet_65mm$decimallatitude)
utms = project(xy, "+proj=utm +zone=11 ellps=WGS84")
Sceloporus_VertNet_65mm$northing = utms[,1]/1000
Sceloporus_VertNet_65mm$easting = utms[,2]/1000
library(sp)
Sceloporus_VertNet_65mm.sp = Sceloporus_VertNet_65mm
coordinates(Sceloporus_VertNet_65mm.sp) = c('easting', 'northing')
Sceloporus_VertNet_65mm.sp$SVL_mm.dev = Sceloporus_VertNet_65mm.sp$SVL_mm - mean(Sceloporus_VertNet_65mm.sp$SVL_mm)
bubble(Sceloporus_VertNet_65mm.sp, zcol = 'SVL_mm.dev', scales = list(draw = T))

library(geoR)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data = Sceloporus_VertNet_65mm$SVL_mm)
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
par(mfrow=c(1,1))
plot(v1)

library(ncf)
ncf.scor <- spline.correlog(Sceloporus_VertNet_65mm$easting, Sceloporus_VertNet_65mm$northing, Sceloporus_VertNet_65mm$SVL_mm,
                            resamp=10, quiet = FALSE)
## 1  of  10 
2  of  10 
3  of  10 
4  of  10 
5  of  10 
6  of  10 
7  of  10 
8  of  10 
9  of  10 
10  of  10 
plot(ncf.scor)

mod.nocorr = lm(SVL_mm ~ bio_6, data = Sceloporus_VertNet_65mm)
v1 <- variog(coords = Sceloporus_VertNet_65mm[,c('easting', 'northing')], data =resid(mod.nocorr))
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
plot(v1)

# jitter GPS points
set.seed(9999)
Sceloporus_VertNet_65mm$j.easting = jitter(Sceloporus_VertNet_65mm$easting)
Sceloporus_VertNet_65mm$j.northing = jitter(Sceloporus_VertNet_65mm$northing)
## corExp
SVL_mm.corExp = gls(SVL_mm ~ bio_6, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corExp(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corExp_4a = SVL_mm.corExp
summary(SVL_mm.corExp_4a)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3728.855 3750.627 -1859.428
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 35.0221945  0.5832795 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.29093 1.2020995 61.80098  0.0000
## bio_6       -0.48914 0.1829269 -2.67394  0.0077
## 
##  Correlation: 
##       (Intr)
## bio_6 0.468 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.88889142 -0.75134844 -0.08037974  0.53208840  4.53643131 
## 
## Residual standard error: 7.338495 
## Degrees of freedom: 577 total; 575 residual
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_4a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

## corGaus
SVL_mm.corGaus = gls(SVL_mm ~ bio_6, 
                     data = Sceloporus_VertNet_65mm, 
                     correlation = corGaus(form = ~j.easting + j.northing, 
                                           nugget = TRUE, 
                                           value = c(50, 0.5)))
SVL_mm.corGaus_4b = SVL_mm.corGaus
summary(SVL_mm.corGaus_4b)      
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3733.052 3754.824 -1861.526
## 
## Correlation Structure: Gaussian spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 28.6290649  0.6175309 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.12319 1.0599616 69.93007  0.0000
## bio_6       -0.53011 0.1636516 -3.23927  0.0013
## 
##  Correlation: 
##       (Intr)
## bio_6 0.476 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.91497565 -0.77036119 -0.08412568  0.54836605  4.56258535 
## 
## Residual standard error: 7.307147 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corGaus.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corGaus_4b, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corGaus.var)

# corLin
SVL_mm.corLin = gls(SVL_mm ~bio_6, 
                    data = Sceloporus_VertNet_65mm, 
                    correlation = corLin(form = ~j.easting + j.northing, 
                                         nugget = TRUE, 
                                         value = c(50, 0.5)))
SVL_mm.corLin_4c = SVL_mm.corLin
summary(SVL_mm.corLin_4c)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3727.277 3749.049 -1858.638
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.4908766  0.6021569 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.25792 1.0075366 73.70246  0.0000
## bio_6       -0.50057 0.1709142 -2.92878  0.0035
## 
##  Correlation: 
##       (Intr)
## bio_6 0.534 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.92996024 -0.76428663 -0.08470616  0.54689828  4.61822848 
## 
## Residual standard error: 7.208484 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corLin.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corLin_4c, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corLin.var)

# corRatio
SVL_mm.corRatio = gls(SVL_mm ~ bio_6, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corRatio(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corRatio_4d = SVL_mm.corRatio
summary(SVL_mm.corRatio_4d)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3733.367 3755.139 -1861.683
## 
## Correlation Structure: Rational quadratic spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 26.6737461  0.6093649 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.09573 1.2292637 60.27652  0.0000
## bio_6       -0.52290 0.1729929 -3.02267  0.0026
## 
##  Correlation: 
##       (Intr)
## bio_6 0.436 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.89383782 -0.75327970 -0.07532492  0.55633844  4.54478478 
## 
## Residual standard error: 7.34641 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corRatio.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corRatio_4d, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corRatio.var)

# corSpher
SVL_mm.corSpher = gls(SVL_mm ~ bio_6, 
                      data = Sceloporus_VertNet_65mm, 
                      correlation = corSpher(form = ~j.easting + j.northing, 
                                             nugget = TRUE, 
                                             value = c(50, 0.5)))
SVL_mm.corSpher_4e = SVL_mm.corSpher
summary(SVL_mm.corSpher_4e)                            
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3729.347 3751.119 -1859.674
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##     range    nugget 
## 47.940820  0.609589 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.22969 1.0122070 73.33450  0.0000
## bio_6       -0.51405 0.1707524 -3.01051  0.0027
## 
##  Correlation: 
##       (Intr)
## bio_6 0.529 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.9480198 -0.7768292 -0.0899644  0.5539295  4.6365309 
## 
## Residual standard error: 7.177685 
## Degrees of freedom: 577 total; 575 residual
SVL_mm.corSpher.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corSpher_4e, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corSpher.var)

# calculate AIC to select best model
AIC(SVL_mm.corExp_4a, SVL_mm.corGaus_4b, SVL_mm.corLin_4c, SVL_mm.corRatio_4d, SVL_mm.corSpher_4e)
# SVL_mm.corLin
# best support
summary(SVL_mm.corLin_4c)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3727.277 3749.049 -1858.638
## 
## Correlation Structure: Linear spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 33.4908766  0.6021569 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.25792 1.0075366 73.70246  0.0000
## bio_6       -0.50057 0.1709142 -2.92878  0.0035
## 
##  Correlation: 
##       (Intr)
## bio_6 0.534 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.92996024 -0.76428663 -0.08470616  0.54689828  4.61822848 
## 
## Residual standard error: 7.208484 
## Degrees of freedom: 577 total; 575 residual
# also good
summary(SVL_mm.corExp_4a)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3728.855 3750.627 -1859.428
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##      range     nugget 
## 35.0221945  0.5832795 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.29093 1.2020995 61.80098  0.0000
## bio_6       -0.48914 0.1829269 -2.67394  0.0077
## 
##  Correlation: 
##       (Intr)
## bio_6 0.468 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.88889142 -0.75134844 -0.08037974  0.53208840  4.53643131 
## 
## Residual standard error: 7.338495 
## Degrees of freedom: 577 total; 575 residual
summary(SVL_mm.corSpher_4e)
## Generalized least squares fit by REML
##   Model: SVL_mm ~ bio_6 
##   Data: Sceloporus_VertNet_65mm 
##        AIC      BIC    logLik
##   3729.347 3751.119 -1859.674
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~j.easting + j.northing 
##  Parameter estimate(s):
##     range    nugget 
## 47.940820  0.609589 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 74.22969 1.0122070 73.33450  0.0000
## bio_6       -0.51405 0.1707524 -3.01051  0.0027
## 
##  Correlation: 
##       (Intr)
## bio_6 0.529 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.9480198 -0.7768292 -0.0899644  0.5539295  4.6365309 
## 
## Residual standard error: 7.177685 
## Degrees of freedom: 577 total; 575 residual
# same result
reg2 <- lm(Sceloporus_VertNet_65mm$SVL_mm~Sceloporus_VertNet_65mm$bio_6)
summary(reg2)
## 
## Call:
## lm(formula = Sceloporus_VertNet_65mm$SVL_mm ~ Sceloporus_VertNet_65mm$bio_6)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.457  -5.200  -0.272   4.309  34.320 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   72.24936    0.56118 128.746   <2e-16 ***
## Sceloporus_VertNet_65mm$bio_6 -0.83810    0.09544  -8.781   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.961 on 575 degrees of freedom
## Multiple R-squared:  0.1182, Adjusted R-squared:  0.1167 
## F-statistic: 77.11 on 1 and 575 DF,  p-value: < 2.2e-16
#variogram
SVL_mm.corExp.var <- variog(coords = Sceloporus_VertNet_65mm[,c('j.easting', 'j.northing')], data = resid(SVL_mm.corExp_4a, type = "normalized"))
## variog: computing omnidirectional variogram
plot(SVL_mm.corExp.var)

# plot
library(ggplot2)
plot2 = ggplot(Sceloporus_VertNet_65mm, aes(x = bio_6, y = SVL_mm)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "black")
plot2 + theme_classic() + labs(x = "Minimum Temperature of Coldest Month", y = "Snout-vent Length (mm)")

# Scales

# Load libraries
library(readr)
library(ggplot2)
library(nlme)
library(plyr)
library(gplots)
library(gstat)
library(sp)

# SE function
se<-function(x) sqrt(var(x, na.rm=T)/length(x))
Sceloporus_specimens = read_delim("Sceloporus_specimens.txt", 
                                  "\t", escape_double = FALSE, trim_ws = TRUE)
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   SeriesID = col_character(),
##   Genus = col_character(),
##   Species = col_character(),
##   Date = col_character(),
##   Area = col_character(),
##   Site = col_character(),
##   Time = col_time(format = ""),
##   Sex = col_character(),
##   Age10g = col_character(),
##   VerbatimLocation = col_character()
## )
## See spec(...) for full column specifications.
summary(Sceloporus_specimens)
##    SeriesID              LACM           Genus             Species         
##  Length:217         Min.   :188168   Length:217         Length:217        
##  Class :character   1st Qu.:188222   Class :character   Class :character  
##  Mode  :character   Median :188278   Mode  :character   Mode  :character  
##                     Mean   :188278                                        
##                     3rd Qu.:188333                                        
##                     Max.   :188387                                        
##                                                                           
##      Date               Area               Site          
##  Length:217         Length:217         Length:217        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##                                                          
##                                                          
##      Time              Sex               Latitude       Longitude     
##  Length:217        Length:217         Min.   :34.10   Min.   :-122.3  
##  Class1:hms        Class :character   1st Qu.:34.37   1st Qu.:-118.6  
##  Class2:difftime   Mode  :character   Median :37.26   Median :-118.4  
##  Mode  :numeric                       Mean   :36.65   Mean   :-118.7  
##                                       3rd Qu.:37.36   3rd Qu.:-117.7  
##                                       Max.   :41.28   Max.   :-116.9  
##                                                                       
##      PosErr        Elevation        SVL_mm         Age10g         
##  Min.   :2.000   Min.   : 973   Min.   :42.00   Length:217        
##  1st Qu.:2.000   1st Qu.:1239   1st Qu.:66.00   Class :character  
##  Median :2.000   Median :1532   Median :73.00   Mode  :character  
##  Mean   :2.451   Mean   :1655   Mean   :72.55                     
##  3rd Qu.:3.000   3rd Qu.:2216   3rd Qu.:79.00                     
##  Max.   :5.000   Max.   :2456   Max.   :93.00                     
##  NA's   :2                                                        
##      Mass_g     VerbatimLocation    Color_Dorsal   Scales_dorsal  
##  Min.   : 2.9   Length:217         Min.   :18.98   Min.   :36.00  
##  1st Qu.:10.3   Class :character   1st Qu.:33.51   1st Qu.:39.00  
##  Median :14.7   Mode  :character   Median :36.75   Median :41.00  
##  Mean   :14.8                      Mean   :36.80   Mean   :41.42  
##  3rd Qu.:17.7                      3rd Qu.:40.49   3rd Qu.:43.00  
##  Max.   :32.1                      Max.   :48.78   Max.   :48.00  
##                                    NA's   :8       NA's   :76     
##  BIO12_annualPrecip BIO5_maxTwarmestMo BIO6_minTcoldestMo Aridity_index   
##  Min.   : 161.0     Min.   :22.4       Min.   :-8.500     Min.   : 13.08  
##  1st Qu.: 365.0     1st Qu.:25.0       1st Qu.:-6.500     1st Qu.: 47.97  
##  Median : 548.0     Median :28.6       Median :-4.800     Median : 72.21  
##  Mean   : 628.4     Mean   :28.9       Mean   :-4.573     Mean   : 87.11  
##  3rd Qu.: 743.0     3rd Qu.:31.6       3rd Qu.:-2.100     3rd Qu.:116.85  
##  Max.   :1760.0     Max.   :35.7       Max.   : 0.400     Max.   :242.62  
##                                                                           
##  BIO10_meanTwarmestQ BIO11_meanTcoldestQ
##  Min.   :13.20       Min.   :-2.200     
##  1st Qu.:16.00       1st Qu.: 1.200     
##  Median :18.10       Median : 2.300     
##  Mean   :18.57       Mean   : 2.447     
##  3rd Qu.:21.70       3rd Qu.: 3.900     
##  Max.   :25.10       Max.   : 7.500     
## 
# Select lizards > 10 g only
Sceloporus_specimens_Age10g = Sceloporus_specimens[which(Sceloporus_specimens$Age10g=='Adult'),]
# select lizards with scale data only
df.complete = subset(Sceloporus_specimens_Age10g, !(is.na(Scales_dorsal)))

plot(df.complete$Scales_dorsal, df.complete$SVL_mm)
SVL.Scales = lm(df.complete$SVL_mm~df.complete$Scales_dorsal)
abline(SVL.Scales, col = "red")

SVL.Scales.resid = resid(SVL.Scales)
SVL.Scales.predict = predict(SVL.Scales)
length(resid(SVL.Scales))
## [1] 107
length(predict(SVL.Scales))
## [1] 107
par(mfrow=c(1,1))
plot(predict(SVL.Scales), resid(SVL.Scales))
abline(0,0)

# POSITIVE residuals are larger scales, NEGATIVE resuduals are smaller scales
resid.Scales_dorsal = resid(SVL.Scales)
df.complete$resid.Scales_dorsal = resid.Scales_dorsal
# Correlations among scale traits

df.complete$Scales_size = (df.complete$SVL_mm/df.complete$Scales_dorsal)

df.complete$Scale.count = df.complete$Scales_dorsal
df.complete$Scale.size = df.complete$Scales_size

# size vs number of scales
cor.test(df.complete$SVL_mm, df.complete$Scales_dorsal)
## 
##  Pearson's product-moment correlation
## 
## data:  df.complete$SVL_mm and df.complete$Scales_dorsal
## t = 4.5891, df = 105, p-value = 1.238e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2372919 0.5554874
## sample estimates:
##      cor 
## 0.408736
# size vs mean scale size
cor.test(df.complete$SVL_mm, df.complete$Scale.size)
## 
##  Pearson's product-moment correlation
## 
## data:  df.complete$SVL_mm and df.complete$Scale.size
## t = 10.721, df = 105, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6178455 0.8026060
## sample estimates:
##       cor 
## 0.7229049

SIMPLE CORRELATION FIGURES

library(ggplot2)
modelR1 = ggplot(df.complete, aes(x = SVL_mm, y = Scales_dorsal)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "black")
modelR1 + theme_classic() + labs(x = "Snout-Vent Length (mm)", y = "Number of Dorsal Scales")

library(ggplot2)
modelR2 = ggplot(df.complete, aes(x = SVL_mm, y = Scale.size)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "black")
modelR2 + theme_classic() + labs(x = "Snout-Vent Length (mm)", y = "Mean Dorsal Scale Size (mm)")

# Mean Temperature of Warmest Quarter

reg3d2.lmer.BIO10 = lmer(resid.Scales_dorsal ~ BIO10_meanTwarmestQ + (1|BIO10_meanTwarmestQ), REML = F, na.action=na.omit, data=df.complete)
#summary(reg3d2.lmer.BIO10)
anova(reg3d2.lmer.BIO10)
plot(reg3d2.lmer.BIO10)

r.squaredGLMM(reg3d2.lmer.BIO10)
##              R2m       R2c
## [1,] 0.003235912 0.1094919
library(ggplot2)
plot2 = ggplot(df.complete, aes(x = BIO10_meanTwarmestQ, y = SVL_mm)) + 
  geom_point()
plot2 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter (°C)", y = "Snout-Vent Length (mm)")

# Maximum Temperature of Warmest Month

reg3d2.lmer.BIO5 = lmer(resid.Scales_dorsal ~ BIO5_maxTwarmestMo + (1|BIO5_maxTwarmestMo), REML = F, na.action=na.omit, data=df.complete)
#summary(reg3d2.lmer.BIO5)
anova(reg3d2.lmer.BIO5)
plot(reg3d2.lmer.BIO5)

r.squaredGLMM(reg3d2.lmer.BIO5)
##              R2m       R2c
## [1,] 0.003710512 0.1119964
library(ggplot2)
plot2 = ggplot(df.complete, aes(x = BIO5_maxTwarmestMo, y = SVL_mm)) + 
  geom_point()
plot2 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month (°C)", y = "Snout-Vent Length (mm)")

# Aridity

reg3d2.lmer = lmer(resid.Scales_dorsal ~ Aridity_index + (1|Aridity_index), REML = F, na.action=na.omit, data=df.complete)
summary(reg3d2.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: resid.Scales_dorsal ~ Aridity_index + (1 | Aridity_index)
##    Data: df.complete
## 
##      AIC      BIC   logLik deviance df.resid 
##    686.9    697.5   -339.4    678.9      103 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.06675 -0.58632  0.07179  0.56288  2.25282 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  Aridity_index (Intercept)  1.489   1.220   
##  Residual                  32.000   5.657   
## Number of obs: 107, groups:  Aridity_index, 30
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)    2.044325   1.013973 22.133786   2.016   0.0561 .
## Aridity_index -0.024633   0.009427 23.243250  -2.613   0.0155 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Aridity_ndx -0.792
anova(reg3d2.lmer)
r.squaredGLMM(reg3d2.lmer)
##             R2m       R2c
## [1,] 0.07119116 0.1124801
library(ggplot2)
reg3d2 = ggplot(df.complete, aes(x = Aridity_index, y = resid.Scales_dorsal)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "black")
reg3d2 + theme_classic() + labs(x = "Aridity Index", y = "Dorsal Scale Size (resid.)")

# Dorsal color
par(mfrow=c(1,1))
hist(Sceloporus_specimens_Age10g$Color_Dorsal, breaks = 20)

# Summer
# Mean Temp of Warmest Quarter

par(mfrow=c(1,1))
par(mar=c(5,5,2,2))

library(lme4)
library(lmerTest)
reg1d2.lmer = lmer(Color_Dorsal ~ BIO10_meanTwarmestQ + (1|BIO10_meanTwarmestQ), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1d2.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Color_Dorsal ~ BIO10_meanTwarmestQ + (1 | BIO10_meanTwarmestQ)
##    Data: Sceloporus_specimens_Age10g
## 
##      AIC      BIC   logLik deviance df.resid 
##    933.0    945.2   -462.5    925.0      153 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.63890 -0.47611 -0.01983  0.61106  2.09243 
## 
## Random effects:
##  Groups              Name        Variance Std.Dev.
##  BIO10_meanTwarmestQ (Intercept) 10.86    3.296   
##  Residual                        16.46    4.057   
## Number of obs: 157, groups:  BIO10_meanTwarmestQ, 33
## 
## Fixed effects:
##                     Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)          29.0645     3.9639 28.2431   7.332 5.24e-08 ***
## BIO10_meanTwarmestQ   0.4148     0.2071 28.6045   2.003   0.0548 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BIO10_mnTwQ -0.984
anova(reg1d2.lmer)
plot(reg1d2.lmer)

r.squaredGLMM(reg1d2.lmer)
##             R2m       R2c
## [1,] 0.06950942 0.4394815
library(ggplot2)
reg1d2 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO10_meanTwarmestQ, y = Color_Dorsal)) + 
  geom_point()
reg1d2 + theme_classic() + labs(x = "Mean Temperature of Warmest Quarter (°C)", y = "Dorsal Lightness")
## Warning: Removed 7 rows containing missing values (geom_point).

# Max Temp Warmest Month

library(ggplot2)
reg1e2 = ggplot(Sceloporus_specimens_Age10g, aes(x = BIO5_maxTwarmestMo, y = Color_Dorsal)) + 
  geom_point()
reg1e2 + theme_classic() + labs(x = "Maximum Temperature of Warmest Month (°C)", y = "Dorsal Lightness")
## Warning: Removed 7 rows containing missing values (geom_point).

library(lme4)
library(lmerTest)
reg1e2.lmer = lmer(Color_Dorsal ~ BIO5_maxTwarmestMo + (1|BIO5_maxTwarmestMo), REML = F, na.action=na.omit, data=Sceloporus_specimens_Age10g)
summary(reg1e2.lmer)
## Linear mixed model fit by maximum likelihood . t-tests use
##   Satterthwaite's method [lmerModLmerTest]
## Formula: Color_Dorsal ~ BIO5_maxTwarmestMo + (1 | BIO5_maxTwarmestMo)
##    Data: Sceloporus_specimens_Age10g
## 
##      AIC      BIC   logLik deviance df.resid 
##    938.7    950.9   -465.3    930.7      153 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9847 -0.4566 -0.0250  0.5928  2.3138 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  BIO5_maxTwarmestMo (Intercept) 10.42    3.228   
##  Residual                       17.15    4.141   
## Number of obs: 157, groups:  BIO5_maxTwarmestMo, 36
## 
## Fixed effects:
##                    Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)         32.0813     4.8770 31.7622   6.578 2.13e-07 ***
## BIO5_maxTwarmestMo   0.1714     0.1641 32.5807   1.044    0.304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## BIO5_mxTwrM -0.990
anova(reg1e2.lmer)
plot(reg1e2.lmer)

r.squaredGLMM(reg1e2.lmer)
##             R2m       R2c
## [1,] 0.01758458 0.3888856

end